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Abstract 

We simulate the lattice Nambu-Jona-Lasinio (NJL) model in 3+1-dimensions at non-zero baryon 
chemical potential Qu) and zero temperature (T) and treat the results as phenomenologically rel- 
evant for cold, dense quark matter. Measurements of the chiral condensate indicate a crossover 
in the thermodynamic limit, whilst at high chemical potential and zero temperature we observe a 
non-zero diquark condensate and a gap in the fermion dispersion relation, which together provide 
evidence for BCS superfluidity. In particular, the size of gap is found to be approximately 15% the 
value of the vacuum fermion mass and roughly independent of fi in the chirally restored phase. 



* Talk presented at the Workshop on QCD in Extreme Environments, Argonne National Laboratory, IL, 
USA, 29 th June to 3 rd July, 2004. 
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I. INTRODUCTION 



Whilst many of the talks at this workshop have been concerned with QCD at extreme 
temperatures, the "extreme environment" in the work presented here is that of cold, dense 
matter, where QCD is believed to exhibit colour superconductivity via the condensation 
of diquarks analogous to the Cooper pairs of BCS superconductors (for recent reviews on 
colour superconductivity see e.g. Q,!^). The BCS mechanism in a superconductor is a subtle 
one, since the fundamental interaction between electrons is repulsive and the net attractive 
interaction due to phonon exchange persists only at extremely low temperatures. In QCD, 
however, the fundamental interaction between quarks is attractive in the anti-triplet channel, 
such that diquark condensation should be far more robust against thermal fluctuations. In 
recent years, studies of four-Fermi models with QCD instanton motivated interactions have 
suggested that the BCS gap in QCD could be as large as 50-100MeV jjj. This means that 
a colour superconducting phase could be relevant to the physics of compact stars where the 
typical temperature is only 0(lMeV). Analytic studies have shown that the ground-state of 
2 + 1 flavour QCD at asymptotically high density is the colour-flavour-locked phase |4j. Due 
to the persistence of the sign problem in lattice QCD with // ^ 0, however, to investigate 
the nature of QCD at more moderate densities one must to resort to studying model field 
theories. The non-perturbative treatment of such models can potentially provide a robust 
method to study colour superconducting matter at the intermediate densities relevant to 
the physics of compact stars. 

QCD with an SU(2) gauge group is one interesting model which can be studied non- 
perturbatively as it does not suffer from a sign problem jf]]. Lattice studies in recent years 
have shown that two colour QCD exhibits superfluidity via Bose-Einstein condensation in 
the dense phase in a manner analogous to He 4 , due to the bosonic nature of the two-quark 
baryons f fl 3. in „ r de r to observe a BCS-st yl e scenario, however, It see m s that one ,nn 8 t 
resort to studying purely fermionic field theories such as the Nambu-Jona-Lasinio (NJL) 
model 0. This model has no gauge degrees of freedom and all orders of gluon exchange 
are approximated by a point-like four-fermion interaction. The original motivation for the 
formulation of this model was that it observed the same global symmetries as strongly 
interacting matter and its vacuum structure exhibited chiral symmetry breaking in a manner 
directly analogous to BCS superconductivity. Therefore, whilst it cannot teach us anything 
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about the mechanism of (de) confinement or the interaction of charged diquarks, if one 
constricts oneself to studying its global symmetries and the pattern of their breaking, it can 
be an ideal model with which to study colour superconductivity in a relativistic quantum 
field theory. 

This talk contains results from a numerical study of the high //, low T phase of the 3 + 1- 
dimensional N JL model, published in , with the aim of showing that the ground-state of 
this model is that of a conventional BCS superfluid formed in a manner analogous to super- 
fluid He 3 . After describing the formulation of the model on the lattice in Sec. |Hj we present 
measurements of the chiral and diquark condensates and discuss the zero temperature and 
infinite volume limits in Sec. 11111 Finally, we present a direct measurement of the BCS gap 
in Sec. II V I and show that it is consistent with model predictions. 



II. THE LATTICE NJL MODEL 

The action of the lattice NJL model, with the lattice spacing a — > 1, can be written as 1 
S = S ferm + S bos = W V^A^y + - 2 £ tr$t* 4 , (1) 

xy 9 x 

where the bispinor \I/ defined on lattice sites x is written in terms of independent isospinors 
via \I/ tr = (x, X tT ) an d the auxiliary scalar and pseudoscalar fields defined on dual lattice 
sites x are introduced via the 2x2 matrix $ = a + m.r. The kinetic operator A written in 
the Nambu-Gor'kov basis is given by 

^ -M tr jr 2 ) 

where j and j are U(1)b symmetry breaking terms, which fix the direction of symmetry 
breaking to allow the measurement of a diquark condensate on a finite volume lattice in 
analogy with the bare quark mass in the measurement of the chiral condensate. In practice, 
we use j and j real, positive and equal, and what we refer to as j from this point on is 
actually the sum of these two terms. Also, simulations are performed in the "partially 
quenched" approximation, with j being zero during the generation of our background field 

1 Here, the bosonic part is corrected from that of the equivalent expression in ^(J by a factor of 2, which 
appeared therein due to a typographical error. 
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configurations, and being made non-zero only during the measurement of the observables. 
This makes the simulations far less computationally expensive, and it should be noted that 
in simulations of the model in 2 + 1-dimensions, there was no discernible difference between 
measurements made using this method and those made via full simulation 
is the standard fermion kinetic operator 



llj. Finally, M 



M pq = \§pq 



v=l,2,3 



+ Jq S *u E W 5 K 9 + ^{x)H{x).r^) , (3) 

(x,x) 

with the parameters being the bare fermion mass m, coupling g 2 and baryon chemical 
potential p. The symbol (x,x) denotes the set of 16 dual sites adjacent to x. The Pauli 
matrices acting on isospin indices p, q are normalised so that trrjTj = 25y. The phase factors 
rj u (x) = (— 1) 2:0-1 ^"-i and e(x) = (—i s ) Xo+Xl+X2+X3 ensure that fermions with the correct 
Lorentz covariance properties emerge in the continuum limit, and that in the limit m —>■ 
the action (0) has a global SU(2)^(S)SU(2)r invariance. In addition, in the limit that j — > 0, 
the action has a U(1)b invariance under baryon number rotations. 

Dimensional analysis shows that the coupling g 2 has mass dimension —2, which reflects 
the fact that in 3 + 1-dimensions, four- fermion models have no non-trivial continuum limit. 
We deal with this by fitting our model to low energy, vacuum phenomenology and extracting 
the relevant bare parameters In particular, we calculate dimensionless ratios of the 

constituent fermion mass and the mass and decay rate of the pion. Setting these to 400MeV, 
138MeV and 93MeV respectively fixes our bare parameters as am = 0.006 and a 2 /g 2 = 0.495. 
For details, see 



III. PHASE STRUCTURE AND THE ZERO TEMPERATURE LIMIT 

Self-consistent treatment of the NJL model shows that for sufficiently strong coupling, 
the approximate SU(2)i®SU(2)^ chiral symmetry is spontaneously broken in the vacuum to 
SU(2)v, leading to a dynamically generated quark mass £ ^> m, and 3 degenerate pseudo- 
Goldstone modes identified with the pions. In the presence of a baryon chemical potential 
fi the symmetry is approximately restored as /x is increased through some onset scale 
p ~ X, with the order of the transition being sensitive to the parameters employed 12 1. 
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We determine the nature of this transition in the regime set by our phenomenological 
parameter choice by studying the order parameter of chiral symmetry breaking, the chiral 
condensate (xx)i defined by 

where Z = J d^d^d^e -5 is the partition function associated with (JI}. (xx) is measured 
for a range of fi over 500 equilibrated hybrid Monte Carlo trajectories of length 1.0 on 
lattice volumes V = L z s x L t = 12 4 , 16 4 and 20 4 with the trace taken via 5 stochastic 
estimations made on every other configuration. The data were then extrapolated linearly in 
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FIG. 1: Chiral and diquark condensates as functions of chemical potential. 

V^ 1 to the infinite volume limit. The data are presented as the circular points in Fig. As 
expected, we see that chiral symmetry is broken at zero chemical potential, as signified by a 
large chiral condensate, and is approximately restored after the condensate passes through a 
broad crossover and becomes approximately zero in the high \i phase. We can also calculate 
(xx) analytically via a self-consistency equation (the gap equation) in the limit that the 
number of fermion species (which we rather arbitrarily refer to as the number of colours N c ) 
is taken to infinity This is plotted as the solid curve and can be seen to qualitatively agree 
with the lattice data, in which the 0(1/N C ) corrections are ~ 15%. 

In order to explore the possibility of a U(l)fi-violating BCS phase at high /i we study the 
relevant order parameter, the diquark condensate (qq+), which in analogy with the chiral 
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condensate is denned by 



1 d\nZ 1 



AV 



V dj 



tr 



T2 



-l 



(5) 



r-2 



Here, the positive subscript highlights the fact that (qq+) is the sum of contributions from 
both XX an d XX states. This expectation value is measured for ten values of j between 
0.1 and 1.0 on the same volumes and with the same parameters with which we measure 
(xx)- Unlike (xx)i however, the data are found not to scale with inverse volume, but 
instead appear to scale linearly with inverse temporal extent. Accordingly, we extrapolate 
our data to L^ 1 — > corresponding to the zero temperature limit. Some of these data are 
presented in Fig. El In contrast with the measurement of (xx)i i n which the bare quark mass 
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FIG. 2: Diquark condensate as a function of j in the zero temperature limit for various \i. 



is left non-zero as the chiral symmetry of QCD is believed to only ever be approximately 
conserved, we are required to extrapolate the data for (qq+) to the j —*■ limit, in which 
the phenomenologically exact U(l) symmetry related to the conservation of baryon number 
is restored. By empirically fitting a quadratic curve through the /i = data, one finds that 
as j goes to zero so does (qq+), implying that there is no condensation in the vacuum as 
expected. For values of \i in the chirally restored phase, however, the picture is not quite so 
clear. It appears that the data correspond to a non-zero diquark condensate, but one can 
only fit a curve though a subset of the data such that one is required to disregard data with 
j < 0.2. Whilst the resulting curves appear to fit the remaining data well, it is important 
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to be able to justify this omission, especially as the data omitted are those closest to limit 
which one is attempting to reach. 

The most obvious possibility is that the condensate is suppressed due to finite volume 
effects, since diquark condensation is associated with the breaking of a global symmetry; as 
the symmetry breaking source is reduced to zero, the correlation length £ of the fluctuations 
of the order parameter should diverge and become comparable to the size of the lattice. 
Figure|Hl however, shows some results from an extensive finite volume study, which illustrates 
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FIG. 3: Diquark condensate vs. j at /i = 0.8 on Lt = 12 lattices with various spacial volumes. 



that (qq+) shows little or no change for any value of j, even when the spatial volume is altered 
by a factor of 2.5 3 ~ 16. It seems impossible, then, that this is the source of the suppression 
of the condensate at low j. Another possibility is that this suppression is due to having 
poor control over the extrapolation to zero temperature. Figure H] shows the results of a 
study of the model at fi — 0.8 and with various inverse temporal extents, corresponding to 
various non-zero temperatures, as well as the curves fitted to j G [0.3, 1.0] of our // = 0.0 
and /i = 0.8 data after the T — > extrapolation. At the highest temperature studied, for 
which L t = 4, if one performs a quadratic extrapolation through the data one observes that 
(qq+) = in the j '• — > limit and that the curve closely resembles the \x = 0.0 curve at 
zero temperature. This suggests that although we are within the chirally restored phase, 
the temperature with L t = 4 is above the critical temperature for superfluidity for all j. As 
the temperature is decreased (i.e. as the L t is increased to 8) the value of (qq+) at j = 1.0 
immediately approaches the curve fitted to our zero temperature results, whilst for lower 
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FIG. 4: Diquark condensate as a function of j at fi = 0.8 at various temperatures. 

values of j this only occurs at even lower temperatures. This can be understood if one 
considers the fact that the effect of j is to make the condensate robust not only to finite 
volume effects, but also to other perturbations such as thermal fluctuations. For a given 
value of j, therefore, there should be some pseudo-critical temperature T c (j) below which 
the data can be said to be in the superfluid phase, and which in the thermodynamic limit 
should increase monotonically with j. This idea is well illustrated by the data measured 
on the 30 3 x 12 lattice, shown in Fig. An attempted fit to a single quadratic through 
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FIG. 5: Diquark condensate as a function of j at fi = 0.8 on a 30 3 x 12 lattice. 

all values of j is of very poor quality, whilst by choosing to split the data at some suitable 
point and fitting the two regions separately, two very reasonable fits can be obtained. In 
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particular, a fit to the data with j < 0.5 is consistent with = suggesting that 

for these values of j, L t = 12 is above the pseudo-critical temperature. The fit to the 
data with j > 0.6, however, is consistent with (qq+)\j_^Q ^ and agrees almost exactly 
with the fit to data from the 12 4 , 16 4 and 20 4 lattices extrapolated to T — > and then 
j G [0.3, 1.0] — > 0, suggesting that this curve represents the true zero temperature limit 
and justifies our discarding of data with j < 0.2. Whilst a linear T — > extrapolation is 
sufficient to reach this curve for 0.3 < j < 0.5, the condensate at j < 0.3 must be suppressed 
too much for such an extrapolation to be sufficient, i.e. the temperature of the lattice on 
our three lattice volumes must be too high compared with T c (j < 0.3). 

Finally, now that we can trust our T — > and j — > extrapolations, let us look back 
at Fig. m where the diquark condensate is plotted with the chiral condensate as a function 
of \x. Although there is clearly a transition from a phase with no diquark condensation to 
one in which the diquark condensate has a magnitude similar to that of the vacuum chiral 
condensate, this transition is far less pronounced than in the chiral case. (qq+) increases 
approximately as /i 2 , but eventually saturates as \i approaches 1.0 and even decreases past 
fi ~ 1.1. This behaviour is directly related to the geometry of the Fermi surface for a system 
defined on a cubic lattice, the area of which we can calculate in the large- N c free fermion 
limit and is plotted as the dashed curve. In the continuum, therefore, (qq+) should continue 
to rise as /i 2 . 



IV. MEASUREMENT OF THE GAP 



In the previous section, we presented evidence for superfluidity in the form of a non-zero 
diquark condensate at high chemical potential. Here we present more direct evidence for 
a superfluid phase by mapping out the fermion dispersion relation (i.e. the energy E as a 
function of momentum k) and observing a BCS energy gap A about the Fermi momentum. 
One advantage of this order parameter over the diquark condensate is that it is directly 
related to a macroscopic property of the superfluid, the critical temperature T c Also, 
being a global order parameter, in principle it would be possible to measure in a gauge 
theory such as QCD, where according to Elitzur's theorem, one cannot write down a local 
order parameter such as (qq + ) in a gauge invariant way [l^ . 
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In order to extract the dispersion relation we define the time-slice propagator 

g(k;t) = Y^g(0,0;x,t)e- fk - g } (6) 

X 

where 



G(x;y) = A X y 



' A N 

■"■xy ly xy 

AT A 

\ xy -"-xy 



(?) 



is the Gor'kov propagator. As in the original BCS theory 15], the fermionic degrees of 
freedom can be viewed as quasi-particles with energy E relative to the system's Fermi 
energy Ep. In the limit that j — ► 0, the propagation of these quasi-particles is described by 
the "normal" (q(0)q(x)) and (q(0)q(x)) parts of (jSJ, i.e. those that are off-diagonal in the 
Nambu-Gor'kov space and related to M~ l . If the Fermi surface is unstable with respect to 
a BCS ground-state, the quasi-particles nearest to Ep undergo particle-hole mixing and a 
gap appears in the energy spectrum. The propagation of these mixed states is generated by 
the diagonal, or "anomalous" (q(0)q(x)) and (q(0)q(x)) parts of (JHJ). By a combination of 
symmetry constraints and empirical observations we note that the complex matrix Q(k,t) 
contains only two independent parts, one in N xy and one in A xy , which from hereon shall 
be referred to simply as the normal and anomalous propagators and written as N(k, t) and 
A(k,t) respectively. 

We measure N(k, t) and A(k, t) at /i — 0.8 on lattices with L x x L y x L z = 96 x 12 x 12 
and L t = 16, 20 and 24 using standard lattice techniques and extrapolate linearly in Lj l 
to zero temperature. This choice means that by choosing k = (k, 0, 0) in (J0J) we can study 
25 independent momentum modes in the x direction between and 7r/2a. We may then 
extract the energy by fitting them to 

N(k,t) = Ae~ Et + Be~ E ^-^ if t = odd 

. (8) 
N(k, t) = if t = even 

and 

A(k,t) = C(e- Et - e - E ( Lt -V) if t = even 
V 1 V ' (9) 

A(k,t) =0 if t = odd, 

where A, B and C are kept as free parameters, as is the energy E, which as expected is found 
to be the same from both (jHJ and (JUJ). These parameters are then, in turn, extrapolated 
to j — > 0. Quadratic polynomial curves are fitted to the coefficients A(k), B(k) and C(k), 
whilst the energy E(k) is fitted with a straight line. As with the extrapolation of (qq+) in 
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the previous section, the extrapolations appear to smoothly fit the data except for at low j, 
where the discrepancy we have attributed to non-zero temperature persists. Again, for the 
purpose of the extrapolations, we believe we are justified in ignoring points with j < 0.3. 




FIG. 6: Propagator coefficients A, B and C at /i = 0.8 in the zero temperature and source limits. 

Figure El shows the coefficients of the normal and anomalous propagators, after the ex- 
trapolations to zero temperature and diquark source, plotted as functions of momentum 
k. Since the normal propagator was chosen from the N xy part of Q, the forward moving 
signal proportional to the coefficient A(k) relates to the propagation of holes in the Fermi 
sea, whilst the backward moving signal proportional to B(k) relates to the propagation of 
particle excitations above the Fermi surface. Accordingly, we see that the low momentum 
excitations near the centre of the Fermi sphere are dominated by hole degrees of freedom, 
whilst the high momentum excitations are dominated by particles. To excitations at the 
Fermi surface, the propagation of fermions forward and backward in time should be equal, 
such that the point where A(k) and B(k) cross allows us to define the Fermi momentum hp. 
The coefficient C(k) ~ at low momentum, but becomes non-zero in a broad peak about 
the position of hp. This vanishing of the anomalous propagator A(k, t), even in the limit 
that j ; — >• is a signal of particle-hole mixing in the presence of a BCS gap. 

For more direct evidence, let us look at the left panel of Fig. which shows the dispersion 
relation at fi — 0.8 and again extrapolated to T — ► and j G [0.3, 1.0] — > 0. The solid curve 
shows the dispersion relation for free massless staggered fermions, which has two distinct 
branches, the hole branch where E(k) decreases with k corresponding to excitations below 
Ep and the particle branch where E(k) increases with k corresponding to excitations above 
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FIG. 7: The lattice dispersion relation and typical free fermion dispersion relation at fi = 0.8. In 
the right-hand panel the hole branch is plotted as negative. 

Ep. In contrast to this, the dispersion relation from our lattice data shows no discontinuity 
between the two branches, which is another sign of particle-hole mixing. More importantly, 
at no time does the curve pass through E = as there is a distinct gap between this point 
and the minimum; this is the BCS gap A = 0.053(6). This can be seen in a more familiar 
light if we plot the hole branch as negative, as illustrated in the right-hand panel of Fig. [7| 
This makes the free fermion dispersion relation a smooth continuous curve, as one would 
expect, whilst for the NJL dispersion relation this introduces a discontinuity at kp, which 
gives a clear illustration of the BCS gap. In order to present the value of the g ap as a 
dimensionless ratio, we also measure the fermion mass at fi = from the vacuum dispersion 
relation which gives us 

A 'f = - 8 » =0.15(2). (10) 
EO = 0.0) V ' K ' 

Assuming a fermion mass of 400MeV, this implies that A(u = 0.8) ~ 60MeV, consistent 

nn 

with the analytic predictions of |3J, [16J . 

Finally, in order to investigate the ft dependence of A, we determine dispersion relations 
for a range of chemical potentials in 0.50 < fi < 0.85, this time using data from L t = 16 and 
20 only. Whilst this means that the extrapolation to T = is no longer an overdetermined 
problem and requires one to estimate the error, it is worth noting that the resulting dispersion 
relation at fx = 0.8 agrees with that presented in Fig. [7| for all k. By calculating E(k) 
and extracting the minimum, we are able to plot A(/i) in Fig. |H1 where it is compared 
with the diquark condensate (qq+)- Whilst the condensate rises with fx, the gap shows 
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FIG. 8: The gap A as a function of fi compared with (qq+) 



no evidence of any /i dependence in the chirally restored phase and a least-squares fit to 
A = constant has a \ 2 °f only 0.33 per degree of freedom. In combination with Fig. 
this supports the simple-minded picture in which only quark pairs within a shell about 
Ep of thickness 2A, independent of fi, contribute to diquark condensation, resulting in a 
condensate (qq+) oc A/i 2 . 



V. SUMMARY 

In this talk we have seen evidence, first presented in ^(|, for a BCS superfluid phase 
at high \i and low T in the 3 + Id lattice NJL model in the form of a non-vanishing local 
order parameter and, for the first time in the systematic study of a relativistic field theory, 
a direct observation of a BCS gap. Given that the model shares its global symmetries with 
the real world, we believe that we are justified in treating this as phenomenological evidence 
for the existence of a similar colour superconducting phase in full QCD. 



ACKNOWLEDGEMENTS 



The author would like to thank Don Sinclair and the staff at the Argonne National 
Laboratory for organising the workshop and for their financial support and hospitality. This 
work was performed in collaboration with Simon Hands. 



13 



APPENDIX: THE LATTICE FERMI SURFACE 

The conclusions drawn from Figs. Eland El (i.e. that (qq+) is approximately independent 
of spatial volume and that one would need to use L t ^> 20 to reproduce the curves in 
Fig. El without discarding any data) might suggest that one could successfully reach the zero 
temperature limit quite cheaply by performing simulations with L t ^> L s . For this reason 
it is worth making the cautionary point that performing simulations with fi 7^ on such 
lattices can lead to unexpected results. At low temperature, the Fermi-Dirac distribution 
resembles a step function, with the discontinuity about the Fermi momentum smeared out 
across a region 5k ~ T ~ Lj 1 ; in a finite system, the Brillouin zone is discretised into a 
cubic momentum lattice with lattice spacing 27i/(aL s ). For lattices with L t ^> L s , therefore, 
the smearing of the Fermi surface is too fine to be resolved on the coarse momentum lattice. 
One consequence is that when the chemical potential is increased smoothly, the Fermi-Dirac 
distribution changes only when the Fermi surface crosses a momentum mode. If one were 
to study a transition, therefore, such as the chiral crossover illustrated in Fig. ^ the physics 
of the system would be approximately constant except at points where the surface crosses 
a mode; the transition would be turned into a series of steps. This is nicely illustrated in 
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FIG. 9: Expectation value of the scalar field vs. fi on a 10 3 x 48 lattice. 

Fig. El where the expectation value of the scalar field £ = (cr) is plotted as a function of 
pna 10 3 x 48 lattice. The Solid curve is £ solved via the gap equation in the large-iV c 
limit and can be seen to change significantly only when the solutions of the free-Fermion 
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dispersion relation in the infinite volume limit 



, ]Tsin 2 A; Fj + £ 2 



(A.l) 



coincide with one of the lattice momenta. The first three of these points are denoted by the 
vertical lines. The lattice data agree qualitatively with this curve and the discontinuities 
can be seen clearly. For this reason one cannot make the ratio L t /L s arbitrarily large. 
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